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ABSTRACT 

We analyze the three-year WMAP temperature anisotropy data seeking to confirm the power spec- 
trum and likelihoods published by the WMAP team. We apply five independent implementations of 
four algorithms to the power spectrum estimation and two implementations to the parameter estima- 
tion. Our single most important result is that we broadly confirm the WMAP power spectrum and 
analysis. Still, we do find two small but potentially important discrepancies: On large angular scales 
there is a small power excess in the WMAP spectrum (5-10% at £ < 30) primarily due to likelihood 
approximation issues between 13 < i < 30. On small angular scales there is a systematic difference 
between the V- and W-band spectra (f ew percent at i > 300). Recently, the latter discrepancy was 
explained bv iHuffenberger et alJ 1)20061) in terms of over-subtraction of unresolved point sources. As 
far as the low-£ bias is concerned, most parameters are affected by a few tenths of a sigma. The 
most important effect is seen in n s . For the combination of WMAP, Acbar and BOOMERanG, the 
significance of n s ^ 1 drops from ~ 2.7a to ~ 2.3a when correcting for this bias. We propose a few 
simple improvements to the \ow-£ WMAP likelihood code, and introduce two important extensions to 
the Gibbs sampling method that allows for proper sampling of the low signal-to-noise regime. Finally, 
we make the products from the Gibbs sampling analysis publically available, thereby providing a fast 
and simple route to the exact likelihood without the need of expensive matrix inversions. 
Subject headings: cosmic microwave background — cosmology: observations — methods: numerical 



1. INTRODUCTION 

Very recently, the Wilkinson Microwave Anisotropy 
Probe ( WMAP) team released the results from 
three years o f observations of the cosmic mic r owave 
background (iHinshaw et alJ l2006t iPaee et alJ l2006t 
iSpergel et alll200fi|) . In addi tion to improving upo n the 
successful first-year release l)Bennett et alJ l2003al) . this 
new data set includes the first full-sky polar ization maps 
at millimeter wavelengths (|Page et al]l200fi|) . and is des- 
tined to be of great importance for the CMB community 
for many years to come. 

Given the dominant position of the WMAP experiment 
in our current understanding of cosmology, it is impera- 
tive that all of the results are verified by several external 
groups using independent techniques. In this paper we 
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begin these efforts by re-computing the angular power 
spectrum and likelihoods 13 , arguably the two most im- 
portant products of any CMB experiment. Similar work 
was carried out after the first-year WMAP release in 2003 
by several groups (e.g., Slosar, Seljak and Makarov, 2004 
- the importance of foreground marginalization at large 
scales; O'Dwyer et al. 2004 - the importance of exact like- 
lihood/posterior evaluation at large scales; Fosalba and 
Szapudi 2004 - pixel space versus harmonic space compu- 
tations). These analyses significantly contributed to the 
improvements that w ere implemented for the three-year 
WMAP data release ijHinshaw et al.ll2006j) . 

Our philosophy in the present paper is that of multi- 
ple redundancy and cross-checking: For both the angular 
power spectrum and the likelihoods (or posteriors) we 
apply several different algorithms and in some cases dif- 
ferent implementations of the same algorithm. Addition- 
ally, the data were downloaded, prepared and analyzed 
by five independent groups, further reducing the risk of 
both programming errors, simple misunderstandings and 
algorithmic peculiarities. This increases our confidence 
in the final results. 

One question that will be given particular attention 
is the issue of statistical and numerical procedure. Ex- 
amples are, how should a high signal-to-noise but low-£ 
likelihood be regularized at its small scale cut-off in or- 
der to avoid numerical artifacts? And in what £-range 
does a pseudo-spectrum estimator provide an adequate 
approximation to the full posterior, and where should an 
exact approach (such as Gibbs sampling) be preferred? 

A separate question is the problem of foreground mod- 

13 We use version v2pl of the WMAP likelihood in this paper; 
http://lambda.gsfc.nasa.gov/ 



2 



eling and subtraction. This will be considered in greater 
detail in a future publication, and for now we mainly 
adopt the approaches used by the WMAP team, with 
one or two notable exceptions: A po wer spectrum that 
does not rely on external information ijSaha et al]l2006j) 
is established and used as a cross-check on the results 
that are obtained from template-corrected maps, and we 
compare spectra computed for different frequencies in or- 
der to search for frequency-dependent signatures. 

Another major part of the three-year W MAP release 
are m easurements of the CMB polarization l)Paee et al.l 
2006). This is a very complex data set, and requires a 
high degree of algorithmic sophistication for proper in- 
terpretation. The necessary extensions to the methods 
described in the present paper are currently being de- 
veloped, and the results from the corresponding analyses 
will be reported upon as soon as they are completed. In 
the present paper, we focus on temperature information 
only. 

The rest of the paper is organized as follows. In Section 
12 we briefly list the methods we use, and in Section El we 
describe the data. Then, we focus on the \ow-£ part 
of the power spectrum in Section 01 and the high-^ part 
in Section |3J Cosmological parameters are considered 
in Section before drawing conclusions in Section [7| 
Algorithmic details are deferred to two Appendices. 

2. METHODS 

In this paper we are primarily interested in scientific 
results, and less in algorithms as such. In this section 
we only list each method that we use, and provide more 
details in the Appendices where necessary, or with refer- 
ences to earlier papers. 

Low-i likelihood evaluation — A central component for 
several of our methods is evaluation of the full likelihood 
C in pixel-space which is defined by 

-21og£ = d*CT 1 d + log|C|. (1) 
up to a normalization constant. Here d are the observed 
data in a pixelized form and C is the corresponding co- 
variance matrix. Evaluation of this quantity involves 
computing a matrix inverse and determinant, and there- 
fore scales as 0(N^ ix ), N p [ x being the number of pixels 
in the data set. Consequently, likelihood evaluations are 
only feasible at low resolutions, and the data must be de- 
graded prior to analysis. Full details on this operation are 
given in Appendix Regularization of the covariance 
matrix with respect to the Nyquist frequency is an issue 
worth some consideration and we do this by wide-beam 
smoothing and noise addition. This is to be contrasted 
to the approach used by the WMAP likelihood code that 
goes b eyond the Nyquist fr equency to achieve a similar 
effect (Hinsha w et al J 12006). For further comments on 
this issue, see Section 0] and Appendix El 

Maximum-likelihood estimation — The maximum- 
likelihood (ML) power spectrum is simply the one that 
maximizes the likelihood as defined by Equation ^ This 
may be found by any non-linear search algorithm, and 
for the present paper we have used both a quasi-Newton 
and a Sequential Quadratic Programming algorithm for 
this task. We obtain identical results with the two meth- 
ods. Previou s analy ses using simil ar techn iq ues include 
iGorskietafl (119961), iBqnd et all 119981). lEfstathioul 
(|2004j) . lSlosar et all l)2004|) and lHinshaw et all (|2006f> . 



Posterior mapping by MCMC — Alternatively, we may 
choose a Bayesian approach and map out the entire 
posterior distribution P(Ci\d) cx £(C£)P(Ce). Here 
P(Ci) is a prior on the power spectrum, which is cho- 
sen to be uniform in the following. Such mapping may 
be done very conveniently with Markov Chain Monte 
Carlo (MCMC) techniques. Some example C MB applica- 
tions o f MCMC techniques ar e de scribed bvlKnox et alJ 
pOOH) . ILewis fc Bridle! p002]) and lEriksen et alJ (|2006ft 7 

Posterior mapping by Gibbs sampling — As mentioned 
above, a direct likelihood evaluation scales as 0(Np ix ), 
and the two above algorithms are therefore limited to 
low resolution maps. Fortunately, it is possible to cir- 
cumvent this problem through a technique called Gibbs 
sampling that allows for sampling from P(Ci, s|d), s be- 
ing the CMB signal, through the conditional distribu- 
tions P(Ci |s, d) and P(s|CV,d). The conditional nature 
of the algorithm avoids inversion of large S + N matrices 
that are dense in all bases. It thereby achieves the same 
goal as the above MCMC technique with d rastically re- 
duced demand on cqmput ational resources ^Jewell et alJ 
l200llWandelt et alJl200l lEriksen et al.ll2004|) . 

MASTER — Approaching the power spectrum estima- 
tion problem from a fundamentally different angle, the 
pseudo-CV methods (e.g., MASTER; Hivon et al. 2002) 
simply 1) compute the spherical harmonics expansion 
with partial sky data; 2) form a raw quadratic estimate 
of the power spectrum; 3) correct for noise; and 4) finally 
decouple the mode correlations by means of an analyt- 
ically computable coupling kernel. Such methods have 
proved to be very useful tools for CMB analysis, pri- 
marily due to their low computational cost, which again 
allows for massive Monte Carlo simulations. 

MASTER with cross-spectra — For multi-channel experi- 
ments, a straightforward and very useful extension to the 
original MASTER algorithm is simply to include only 
cross-correlatio ns between channels in the power spec- 
trum estimate l|Hinshaw et all 12003). While an auto- 
correlation implementation is quite sensitive to assump- 
tions in the noise model, such assumptions only affect 
the error bars in a cross-spectrum and not the mean. 

MASTER with internal foreground cleaning — For multi- 
frequency experiments with multiple channels per fre- 
quency it is possible to form a set of foreground-cleaned 
maps using different channels within each frequency. The 
noise contributions are thus still uncorrelated among 
several pairs of cleaned maps, and the cross-correlation 
MASTER implementation may be applied as described 
above even to these maps. Su ch a method was im- 
plemented by ( Sa ha et alJ l200 6[) using t he foreground 
cleaning method of iTegmark et all l)2003fl . We use this 
method in the present paper for monitoring residual fore- 
grounds in the spectra computed with template-cleaned 
maps, and will refer to the method as "MASTER with 
internal cleaning" (MASTERint for short; MASTERext 
refers to foreground correction by external templates). 

3. DATA 

The WMAP temperature data l|Hinshaw et al] 1*2006) 
are provided in the form of ten sky maps observed at 
five frequencies between 23 and 94 GHz. These maps 
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are pixelized at HEALPix resolution parameter -/V s id c = 
512 resulting in about 3 million pixels per map. In this 
paper we mainly consider the V-band (61 GHz) and W- 
band (94 GHz) channels since these are the ones with 
the lowest foreground contamination levels. However, 
the full data set is used by MASTERint. 

We take into account the (assumed circularly sym- 
metric) beam profile of each channel independently, and 
we adopt the Kp2 sky cut as our base mask. This 
excludes 15.3% of the sky including all resolved point 
sources. However, for the low-resolution analyses we 
consider downgraded and/or expanded versions of this 
mask. The details of the degradation procedure will be 
discussed later. 

The noise is modeled as uncorrelated, non-uniform and 
Gaussian with an RMS given by <Jo,i/ \/-^obs (p) • Here 
(To,i is the noise per observation for channel i, and Aobs(p) 
is the number of observations in pixel p. However, this 
approximation is not adequate for all channels, and in 
particular the W-band must be treated with special care 
because of significant noise correlations. 

We consider two different foreground correction pro- 
cedures. First, we simply use the template-corrected 
maps as provided on the LAMBDA w ebsite, from which 
a free- free template llFinkbeinerir2003|) . a dust template 
ijFinkbeiner et al.lll999|) and the K-Ka WMAP sky map 
have been fitted and subtracted. Second, to cross-check 
these results we include a MASTER implemen tation that 
does not rely on any external information ijSaha et alJ 
2006). 

Finally, the power spectrum contributions from unre- 
solved point sources are subtracted band-by-band as a 
po st-processing step acc ording to the model described 
bv lHinshaw et aU l|2006fl . 

4. LARGE-SCALE ANALYSIS 

In this section we present the results from our large- 
scale analyses, broadly defined by t < 50. The high- 
resolution results are shown in Section [5j 

4.1. Pre-processing and data selection 

Since pixel-space likelihood evaluations scale as 
0(ATp ix ), we require low-resolution data for both the 
maximum- likelihood and MCMC analyses. For these 
methods, we therefore consider downgraded maps at 
= 16, each having 3072 pixels on the full sky. 
For Gibbs sampling, angular resolution does not carry 
a similar resolution-dependent computational cost, and 
the analysis is made with full-resolution data. MAS- 
TER results are ignored until the next section, since they 
are sub-optimal and difficult to compare with other es- 
timates at large angular scales. However, we note that 
we have reproduced the MASTER spectra presented by 
the WMAP team exactly with MASTERext, and very 
closely with MASTERint. 

The full-resolution sky maps were downgraded as fol- 
lows. (See Appendix [A] for a thorough discussion of 
this procedure.) First, the full resolution maps were de- 
convolved by the instrumental beam and N S id e = 512 
pixel window, and then re-convolved by a Gaussian 9° 
FWHM beam and N s id e = 16 pixel window. Finally, uni- 
form Gaussian noise with a standard deviation of 1 /xK 
was added to each pixel. This ensures that the map is 
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Fig. 1. — Top panel: Comparison of low-£ power spectra com- 
puted from the three-year WMAP data using different techniques 
and data combinations. Bottom panel: Differences between the 
spectra computed in this paper and the (appropriately binned) 
WMAP spectrum. See text for full details. 



properly bandwidth limited at ^ max = 3 • A^dc = 48 and 
noise dominated beyond £ = 40. 

The full-resolution Kp2 mask was downgraded in two 
different ways. In the first case, we set excluded pixels 
to zero and included pixels to unity. Then we downgrade 
the mask using the hierarchical structure of HEALPix, 
and set each low-resolution pixel to the average of all 
sub-pixels. If this average is larger than 0.5, the low- 
resolution pixel is accepted. This is the method used 
in the WMAP likelihood code, and will be denoted by 
"Kp2" in the following. (In order to compare directly 
to the WMAP results, we first downgrade to -/V s id = 8, 
and then upgrade to A S jd c = 16, ensuring identical sky 
masks.) In the second case, we smooth the original mask 
by a Gaussian beam of 9° FWHM, project the smoothed 
field onto the low-resolution grid, and reject all pixels 
with a value smaller than 0.991, roughly expanding the 
original mask by one FWHM in all directions. This mask 
will be denoted "Kp2e" in the following - "e" for ex- 
tended - and excludes about 7% more pixels than the 
Kp2 mask. 

We consider three maps in these analyses, namely the 
de-biased Internal Linear Combination (ILC) map and 
the template-corrected V- and W-band maps individu- 
ally. Instrumental noise is negligible at these scales, and 
co-addition is neither required nor desirable. The non- 
zero frequency range is instead used for foreground mon- 
itoring. 

In order to further minimize the impact of foregrounds, 
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Fig. 2. — Power spectrum dependence on sky cut and data 
set. Colored lines show the maximum likelihood spectra (computed 
with a pixel-based likelihood) for different sky cuts (blue - directly 
downgraded Kp2; red — expanded Kp2) and data sets (dashed - 
ILC; solid - V-band). Note that the expanded Kp2 mask spec- 
tra lies systematically below the directly downgraded Kp2 mask 
spectra. 



we adopt the approach of the WMAP team, and con- 
struct a foreground template as the difference between 
the raw V-band and the ILC map. This template is then 
downgraded in a similar manner as the data maps, but 
without noise addition. A corresponding term is added 
to the covariance matrix with a large pre-factor in order 
to remove any sensitivity to fluctuations with the same 
spatial pattern as the template (e.g., Bond et al. 1998 
and Appendix El . 

4.2. Specification of algorithms 

Three different algorithms are used for the low-^ analy- 
sis, namely Gibbs sampling, Maximum Likelihood (ML) 
estimation with low-resolution data, and MCMC with 
low-resolution data. In the first case, two different im- 
plementations are used (Commander and MAGIC; Erik- 
sen et al. 2004) using different priors. In both cases, the 
power s pectra were refin ed using the Blackwell-Rao es- 
timator l)Chu et al.ll200. l Tfl . For Commander, we use uni- 
form priors, and report the modes of the posterior dis- 
tributions as the power spectrum estimates. This may 
be directly compared to the WMAP spectrum which is 
claimed to be a maximum-likelihood estimate at low ts. 
For MAGIC, we use the Jeffrey's ignorance prior, and 
report the means as our power spectrum estimates. In 
either case, we adopt asymmetric 68% confidence regions 
as our uncertainties. 

The ML estimates are established by means of two dif- 
ferent non-linear search algorithms from a commercial 
library, one quasi-Newton method and one sequential 
quadratic programming method. The free parameters 
are the binned CiS (adopting the binning scheme used 
by the WMAP team) up to I = 30, for a total of 9 
free variables. From I — 31 the spectrum is fixed at the 
WMAP spectrum values. The covariance matrix is de- 
fined as detailed in Appendix and consists of a sum 
of a CMB signal term, a noise term and a foreground 
(monopole, dipole and ILC-V difference map) term. 

The MCMC analysis is identical to the ML analysis 



in terms of the likelihood evaluation method, binning 
scheme, multipole range etc., but uses a Metropolis- 
Hastings algorithm to probe the full posterior distribu- 
tion. We adopt a uniform prior even in this case. Eight 
independent chains are run for 100 000 steps each, en- 
suring excellent convergence properties. One single step 
costs approximately 1 CPU second for a 2000 pixel data 
set. 

4.3. Results 

In Figure ^ we show a comparison of four different 
\ow-£ power spectrum estimates: The WMAP spectrum 
(black), the Commander posterior mode uniform prior 
V-band spectrum (red), the MAGIC posterior mean Jef- 
frey's prior V-band spectrum (orange), and the ML V- 
band spectrum computed with the Kp2e mask (blue). In 
the bottom panel, we show the differences taken between 
these spectra and the WMAP spectrum. Error bars are 
only shown in the top panel: All analyses study the same 
data set, and cosmic variance uncertainties can therefore 
not be used to compare the relative agreement between 
the various solutions. 

First we note that the individual variations between 
all four spectra are about 50-100 //K 2 , corresponding to 
about 5-10% of the absolute spectrum amplitude. This 
is much less than the cosmic variance, and the analy- 
sis uncertainty per mode is thus not a dominant effect. 
However, it is troublesome that this difference appears 
to be systematic: All three colored spectra lie below the 
WMAP spectrum up to I « 30 — 40. Since the cumula- 
tive effect over an entire multipole decade of even a small 
bias may potentially affect cosmological parameters, it is 
important to determine the origin of this discrepancy. 

Foreground uncertainties are always a concern at large 
angular scales for CMB experiments and a most natural 
first candidate to consider. First, we note that the four 
spectra shown in Figure ^ are computed from slightly 
different data sets: The WMAP spectrum is based on the 
ILC map with a directly downgraded (i.e., not expanded) 
Kp2 mask for I < 12 and on the template-corrected V- 
and W-bands with the high- resolution Kp2 cut at £ > 12; 
the Gibbs spectra are based on the template-corrected V- 
bands and the full-resolution Kp2 cut at all scales; the 
ML estimate is computed from the degraded template- 
corrected V-band, cut with the expanded Kp2 mask. 

To assess the impact of residual foregrounds in the 
\ow-£ power spectra, we show in Figure |2] the ML power 
spectra computed from the ILC and template-corrected 
V-band maps, with both the original Kp2 and extended 
Kp2e sky cuts. A noticable trend is clearly visible in 
this range, in that there is a significant loss of power 
between the Kp2 and Kp2c cut for both the ILC and 
V-band maps. This implies that there is considerable 
additional power in the 7% near-galactic pixels in the 
Kp2 sky cut relative to the Kp2e cut for both maps, and 
both maps are most likely contaminated at some level 
near the mask boundaries. On the other hand, after 
expanding the mask the two spectra are quite similar, 
possibly indicating that one is fairly safe after expanding 
the mask, and that the high latitude residuals are small 
compared to the CMB fluctuations. 

While the mask explanation can account for some of 
the discrepancies seen in Figure ^ it only does so at 
I < 12 where a pixel-based likelihood evaluation is used 
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Fig. 3. — Posterior distributions (and the WMAP likelihood) computed by different methods and data sets. The colored histograms 
show the results from MCMC runs with low-^ likelihood computed in pixel space for different masks, and the black curve shows the 
Blackwell-Rao estimator computed with V-band/Kp2 Gibbs sampl es. The vertical dashed lines show the angular power spectrum for the 
ACDM model that best fits the WMAP data (Hinshaw ct al. 2006). The horizontal dashed lines shows the 95% confidence region relative 
to the Gibbs-based estimate. Note that while the MCMC shows marginalized distributions, the other curves indicate slices through the 
multi-dimensional distributions with other Ci's fixed at the best-fit spectrum value. 



by the WMAP team. For I > 12 a MASTER-based 
likelihood is used which is based on the full-resolution 
template-corrected data. The small discrepancies seen 
in the range between 12 < I < 30 is therefore due to dif- 
ferences in statistical treatment, and not data selection. 

A main question is therefore the following: For what 
Grange does the MASTER-based likelihood approxima- 
tion provide an acceptable fit to the exact likelihood? 
From Figure it seems clear that £ e xact < 12 is not suffi- 
cient, while ^exact < 50 is quite likely more than enough. 
It is difficult to answer this question based on a power 
spectrum plot alone, and we will therefore return to this 
question in SectionHOwhere we estimate cosmological pa- 
rameters with various likelihoods. 

In order to further illuminate the above issues, we 
show in Figure |31 a set of different likelihood (and pos- 
terior) distributions computed from the WMAP data. 
The vertical dashed lines show the binned best-fit 
ACDM spectrum and the curves show the WMAP like- 
lihoods (black), Gibbs + Blackwell-Rao posterior distri- 
butions (red/orange) and MCMC posterior distributions 
(blue/green), respectively. Two sigma confidence regions 



relative to the Gibbs distributions with uniform priors 
(red curves) are indicated by horizontal dotted lines. 
Many interesting points may be seen in this figure: 

1. Comparing the V-band+Kp2e and ILC+Kp2 
MCMC runs (green versus blue) we see that the 
former is shifted systematically to lower values of 
Cg. This shift is most likely due to a combination 
of extra fluctuation power coming from the region 
excluded by Kp2e but included by Kp2, and ad- 
ditional sampling variance from the extended sky 
cut. 

2. The agreement between the ILC+Kp2 MCMC and 
WMAP likelihood distributions (blue versus black) 
is generally quite reasonable, but small shifts may 
be seen in a number of cases, especially for 6 < 
I < 12. This is likely due to different likelihood 
regularization (see Appendix \K$ : We use a prop- 
erly bandwidth limited map within the Nyquist fre- 
quency £n = 3-/V s id but add a small noise term 
to regularize the covariance matrix; the WMAP 
code increases ^ ma x to 4A r s i ( j e , well beyond the 




We now consider the intermediate- 
parts of the power spectrum. 



and small-scale 
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Fig. 4. — Top panel: Comparison of power spectra computed 
from the three- year WMAP data using different techniques and 
data combinations. Bottom panel: Differences between the spectra 
computed in this paper and the (appropriately binned) WMAP 
spectrum. 



Nyquist frequency, to make the covariance matrix 
non-singular. 

3. The V-band+Kp2e MCMC distributions are gener- 
ally wider than the ILC+Kp2 distributions (green 
versus blue) because of a more conservative mask. 
This is avoided with the full-resolution Gibbs anal- 
ysis (red curve), since it uses unsmoothed data and 
thereby does not need to expand the mask. Conse- 
quently, the Gibbs distributions are about as wide 
as the ILC+Kp2 distributions. 

4. The Jeffrey's ignorance prior puts more probability 
on small values of Ce, but is mostly relevant at very 
low £'s (orange versus red). Because of the very 
large cosmic variance of these multipoles, the effect 
of the Jeffrey's prior on cosmological parameters is 
small. 

5. The quadrupole posterior maximum with Jeffrey's 
and uniform priors are respectively 105/iK 2 and 
WOfiK 2 . Thus, even with Jeffrey's prior the 
WMAP quadrupole amplitude is completely unre- 
markable relative to the best-fit ACDM value, and 
always well inside the 95% confidence region (ver- 
tical dashed versus horizontal dotted lines). The 
most extreme outlier is that of 20 < £ < 24, whose 
value is low at the 99.4% level. 

5. INTERMEDIATE- AND SMALL-SCALE ANALYSIS 



5.1. Algorithms and data selection 

Since we study full-resolution data in this section, 
the direct likelihood evaluation techniques are no longer 
available to us. However, on these scales the MASTER 
algorithms are applicable, and we once again have mul- 
tiple methods available for cross-checking purposes. In 
total four different codes are used, namely Commander 
(Gibbs sampling; uniform prior), MAGIC (Gibbs sam- 
pling; Jeffrey's prior), MASTERext and MASTERint. 

For the Gibbs analyses we consider only the V-band 
data; the Q-band is considered to be too foreground con- 
taminated for reliable analysis, and the W-band exhibits 
strong correlated no ise that s ignificantly biases any auto- 
correlation method ijHinshaw et alJ l2006) . This has been 
verified in our analyses, and we exclude these bands en- 
tirely from the Gibbs analyses. Next, we analyzed the 
data with both Commander and MAGIC, and obtained 
identical results (up to convergence) on small and inter- 
mediate scales where the prior is less important. There- 
fore we show only one set of Gibbs results in the follow- 
ing. 

In the case of the MASTER analyses, we include cross- 
spectra only in the final power spectra, and are thus quite 
safe with respect to correlated noise. For this reason the 
W-band is included in these cases. Further, for the in- 
ternal cleaning analysis two different data combinations 
were considered, namely including either all five bands 
or only the cleanest Q-, V- and W-bands. The difference 
between these two spectra was found to be very small. 
We once again present only one (the five-band) solution 
here, and note that this solution is not strongly depen- 
dent on the low-frequency bands. 

The foreground mask is always the full-resolution Kp2 
that excludes resolved point sources. Contributions to 
the power spectrum from unresolved point sources are 
subtracted band-by-band accor ding to the model de- 
scribed bv iHinshaw eta!] p00?j) . 

5.2. Results 

A summary of our main intermediate- and small-scale 
power spectra are shown in Figure 0] Four spectra are 
plotted here: The WMAP spectrum, the Gibbs V-band 
spectrum, the V+W MASTERext cross-spectrum, and 
the five-band MASTERint cross-spectrum. 

In the top panel we see that the overall agreement is ex- 
cellent considering the very different approaches taken by 
the different methods. It is thus very unlikely that algo- 
rithmic issues play a dominant role in the determination 
of cosm ological parameters as presented bv lSnergel et alJ 
( 2006) for the small-scale temperature anisotropy results. 

In the bottom panel we plot the difference between 
each of the colored spectra and the WMAP spectrum in 
order to look for systematic relative trends. Again, we 
see that the overall agreement is very good, and the only 
possible anomaly is a slight power deficit for the V-band 
Gibbs spectrum at £ fts 300-600. This will be studied 
further shortly. 

The extreme high-£ discrepancies seen in two MAS- 
TER spectra are due to noise weighting differences. For 
the MASTERext analysis, we obtained the individual 
cross-spectra from the WMAP team and verified that 
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Fig. 5. — Power spectrum versus frequency. All spectra 
are computed from the template-corrected maps with the basic 
cross-MASTER, algorithm. Only the frequency contents vary, as 
each spectrum only contains cross-information of Vx V, VxW and 
WxW respectively. 



our spectra are identical to those of the WMAP team. 
However, while the WMAP spectrum uses an elaborate 
and nearly maximum-likelihood weighting scheme for co- 
adding these cross-spectra into one optimal spectrum, we 
use either simple inverse-variance weighting (for MAS- 
TERext) or flat weighting (for MASTERint). Differences 
in the very low signal-to-noise regime are thus not unex- 
pected. 

The excellent agreement seen in Figure 0] is perhaps 
most remarkable in the case of the MASTERint imple- 
mentation: While all other spectra more or less corre- 
spond to different statistical treatment of the same ba- 
sic data set, this particular solution takes a drastically 
different approach, and uses only WMAP data by them- 
selves to correct for foregrounds. Thus, it provides an im- 
portant and impressive cross-check on the WMAP fore- 
ground cleaning approach. 

However, returning to the power deficit issue seen in 
the V-band spectrum, we plot in Figure [3] the MAS- 
TERext spectra for individual frequency combinations, 
including VxV, VxW, or WxW spectra, respectively. 
In this plot we see that the MASTER V-band spectrum 
agrees very well with the Gibbs V-band spectrum, and 
lies systematically below the averaged WMAP spectrum 
between £ ss 300 and 600. Furthermore, we see that the 
W-band spectrum lies systematically above the WMAP 
spectrum in the same range, with a relative bias of up to 
80 \i¥}. The VxW cross-spectrum lies in between, being 
more or less an average of the two. 

To quantify the significance of this difference, we ana- 



lyzed the set of 2500 Monte Carlo simulations that was 
used for the MASTERext analysis. For each simulation 
we first computed the difference between the W- and 
V-band spectra, and then the mean (inverse noise vari- 
ance weighted) difference within some ^-range. Compar- 
ing the observed WMAP data to these simulations, we 
found that the discrepancy is significant at slightly more 
than 3cr for both 250 < t < 600 and 400 < £ < 500, 
corresponding to the largest region of asymmetry and to 
the most discrepant region, respectively. Presumably it 
is possible to boost the significance by tuning the region 
further, but on the other hand, the known ~ 1% beam 
amplitude uncertainty that is taken into account by the 
WMAP likelihood code will reduce the signficance by a 
few tenths of a sigma. Independent of such minor effects, 
it seems that the probability of this being a statistical 
fluke is rather small. 

After t he publication of the presen t paper, a follow-up 
study by iHuffenberger et all l)2006l) revisited the unre- 
solved point source analysis performed by the WMAP 
team. The main result from that work was a best- 
fit unresolved point source spectrum amplitude of A = 
0.011 nK 2 sr (relative to 41 GHz), which is significantly 
lower t han the value of A = 0.017 /iK 2 sr initially re- 
ported ( Hinsh aw et al.ll2006|) and used in the present pa- 
per. Thus, a relatively larger contribution is subtracted 
from the V-band than from the W-band spectrum, and 
in effect, the V-band spectrum has been over-corrected. 
After taking into account this new amplitude, the dis- 
crepancy between the two band spectra was found to be 
significant at ~ 2a using the same test as above. 

While the quoted point source correction can account 
for a substantial amount of the observed discrepancy, 
there is still a small difference present, and this could 
possibly indicate further systematics. In this respect, 
it could be worth considering unmodeled beam asym- 
me tries. As a prel i minar y step in the three-year analy- 
sis, iHinshaw et al.l l)2006j) considered the impact of such 
asymmetries on the individual cross-spectra. After an 
extensive analysis, they concluded that their magnitude 
is less than 1% at £ < 1000, and therefore sufficiently 
small to neglect in further analyses. However, at the rel- 
evant scales, the absolute amplitude of the temperature 
power spectrum is about 2000 /zK 2 , and a beam asym- 
metry bias of 1% therefore corresponds to an expected 
discrepancy of 20 /iK 2 . This corresponds roughly to la 
of the V- versus W-band difference, and, if determined 
appropriate, its correction could bring the overall differ- 
ence well within the statistical errors. Fortunately, this 
issue will be clearer with additional years of WMAP ob- 
servations. 

6. COSMOLOGICAL PARAMETERS 

In the previous two sections we considered the temper- 
ature angular power spectrum as observed by WMAP, 
and our main conclusion from these analyses is that the 
WMAP spectrum is reasonable at all angular scales. 
However, there are small but clearly noticeable biases 
at both large and small scales. In this section, we 
seek to quantify the impact of this bias in terms of 
the cosmological parameters for a minimal six-parameter 
ACDM model. The combined effect of both the low- and 
high-l discrepancies are studied by IHuffenberger et all 
(2006), and the effect on extended cosmological models 
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Fig. 6. — Probability contours for a simple amplitude/tilt model 
(C e (q,n) = qCf d (l/10) n , Cf d being the best-fit ACMB model) 
fitted to the WMAP likelihood (solid lines) and to the Gibbs sam- 
pled posterior (dotted lines). Contours indicate 68 and 95% confi- 
dence regions. Note the shift to higher amplitudes for the WMAP 
results, and also that the Gibbs results are more closely centered 
on (q,n) = (1,0). 



(e.g., models including massive neutri nos and running of 
the sp ectral index) are considered bv iKristiansen et al.l 
(200j§. 

We perform four sets of similar analyses, all pri- 
marily based on the WMAP likelihood code. First, 
we adopt the WMAP likelihood as provided. Second, 
we replace the \ow-£ likelihood (both the pixel-based 
estimator and the \ow-£ MASTER estimator) with a 
Blackwell-Rao (B R) Gibbs-based estimator for £ < 30 
l|Chu et alJl2005|) . Third, we do the same for £ < 12 
alone. Finally, we use the A^idc = 16, ^ max = 30 
pixel-based likelihood for the V-band and Kp2e cut de- 
scribed earlier. For all cases, we analyze two data com- 
binations; the W MAP data alone, and the combination 
of WMAP. Acbar llKuo et al.ll200l and BOOMERanG 
( Mont rov et alJl200Rt iPiacentini et al]l2005t I.Tones et al.l 
2005). Marginalization of SZ was not included. 

Note that we use the WMAP likelihood at high £'s 
in all cases in order to highlight the effects of the low- 
£ likelihood bias. While using a Gibbs sampling based 
estimator even at high Ps could potentially have benefi- 
cial effects in terms of uncertainties, it might confuse the 
low-£ bias analysis by introducing other differences. 

The cosmological parameters corresponding to these 
likelihoods were established through standard Markov 
Chain Monte Carlo analysis. In keeping with our philos- 
ophy of cross-verification, two independent codes were 
used for a few cases . In the first case, CosmoMC 
ijLewis fc Bridld l2002j) was downloaded and appropri- 
ately modified, and in the second, a stand-alone code was 
written from scratch in C-l — h The two codes returned 
identical distributions, and as usual we show only one 
set of results in the following. We also performed similar 
analyses using only temperature-data, imposing a Gaus- 
sian prior on the optical depth of r = 0.10 ± 0.03 to 
simulate the effect of polarization data, but without re- 
lying on the accuracy of these. As expected, we then 
obtained very similar results to those reported here. 

Before reporting the physical parameters, we consider 



a very simple phenomenological model in order to gain 
intuition on what to expect. Specifically, we fit a model 
with a free amplitude and tilt to both the WMAP like- 
lihood and the WMAP+BR £ max = 30 hybrid, 

C e ( q ,n) = qC? d ^y . (2) 

Here Cf d is a fiducial model (chosen to be the best-fit 
ACDM power law model), q is a relative amplitude, and 
n is a tilt parameter. Since we expect the fiducial model 
to be a good fit to the data, we anticipate values of (q, n) 
around (1,0). In this exercise we include only data be- 
tween 2 < £ < 20. The results from these computations 
are shown in Figure [5] where contours indicate 68 and 
95% confidence regions. WMAP results are indicated by 
solid lines and Blackwell-Rao results by dotted lines. 

By far the most striking feature is a ~ 0.5er shift 
in amplitude towards lower values for our revised pos- 
terior, consistent with the \ow-£ power spectrum re- 
sults shown earlier. However, it is interesting to notice 
that the Blackwell-Rao contours are more centered on 
(q,n) — (1,0) than the WMAP contours: This indicates 
that there is less tension between the low-£ and high-i 
parts of the spectrum when using our approach instead 
of the WMAP approach. 

For the physical parameter estimation, we adopted 
a standard six-parameter ACDM model with 
{fl h h 2 ,D, c h 2 ,9,T,n s ,\og(10 10 A s )} as our free pa- 
rameters. The priors were chosen to be uniform for all 
parameters. 

The results from these computations are summarized 
in Table and Figure [7| and a number of interesting 
points may be seen here. First, comparing the WMAP 
likelihood to the £ < 12 WMAP + BR hybrid (second 
and third columns of Table 0) , we see that the small 
change due to mask expansion at £ < 12 seen in Figure 
121 has no impact in terms of cosmological parameters. In 
effect, the cosmic variance on these very large scales hides 
such problems. 

However, comparing the WMAP likelihood to the £ < 
30 WMAP + BR hybrid (second and fourth columns), 
a different picture is seen. In this case, we find shifts 
of up to 0.4cr for both h and n s (right column). Thus, 
the MASTER-based approximation used in this regime 
appears to be inadequate for the precision required by 
the WMAP data. That the differences are indeed due 
to statistical treatment is confirmed by the results for 
the low-£ pixel-based hybrid, which extends the original 
WMAP analysis simply by using direct evaluation up to 
£ = 30 instead of I = 12. 

We also performed a second analysis with the WMAP 
+ BR hybrid, this time making the transition at £ = 50. 
Fortunately, the results from these computations were 
essentially identical to those for £ = 30, and this implies 
that a pixel-based likelihood for N^ B = 16 may still be 
used for the precision of the WMAP data. However, the 
computational expense of this approach is already push- 
ing close to what is reasonable, and currently dominates 
the MCMC cost. This illustrates very well the advan- 
tages of the Gibbs sampling approach, in which identical 
results are obtained at only a fraction of the computa- 
tional cost. 

In Figure we compare the results for the WMAP 
likelihood and WMAP + BR hybrid likelihood from the 
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TABLE 1 

COSMOLOGICAL PARAMETERS 



WMAP + BR WMAP + BR WMAP + Pixel Shift 
Parameter WMAP (^ max = 12) (£ max = 30) (£ max = 30) (in ctwmap) 



WMAP data only 



Q b h 2 


0.0222 ±0.0007 


0.0222 ± 0.0007 


0.0224 ± 0.0007 


0.0224 ± 0.0007 


-0.3 




0.241 ±0.036 


0.241 ± 0.0037 


0.230 ± 0.033 


0.234 ±0.035 


0.3 


log(10 10 A s ) 


3.019 ±0.067 


3.017 ±0.068 


3.014 ±0.068 


3.013 ±0.068 


0.1 


h 


0.731 ±0.033 


0.731 ±0.032 


0.743 ± 0.033 


0.739 ±0.033 


-0.4 


n a 


0.954 ±0.016 


0.955 ±0.016 


0.961 ± 0.016 


0.960 ±0.017 


-0.4 


T 


0.090 ± 0.030 


0.088 ± 0.030 


0.090 ± 0.030 


0.088 ±0.030 


0.0 






WMAP + Acbar 


+ BOOMERanG 






n h h 2 


0.0225 ± 0.0007 


0.0225 ± 0.0007 


0.0227 ± 0.0007 


0.0227 ± 0.0007 


-0.3 




0.239 ±0.031 


0.240 ± 0.032 


0.229 ±0.030 


0.233 ±0.031 


0.3 


log(10 10 A s ) 


3.030 ± 0.064 


3.028 ± 0.065 


3.025 ± 0.066 


3.023 ±0.0065 


0.1 


h 


0.737 ±0.029 


0.736 ± 0.031 


0.749 ±0.031 


0.744 ±0.031 


-0.4 


n B 


0.958 ±0.016 


0.959 ±0.016 


0.965 ±0.016 


0.964 ±0.016 


-0.4 


T 


0.091 ±0.030 


0.090 ± 0.030 


0.091 ± 0.030 


0.089 ± 0.030 


0.0 



Note. — Comparison of marginalized parameter results obtained from the WMAP likelihood 
(second column), the WMAP + Blackwell-Rao hybrid (third and fourth columns), and from 
the WMAP + Af s idc = 16 pixel-based likelihood (fifth column) . The latter three were based on 
the template-corrected V-band data at low ^'s. The relative shift between the WMAP and the 
^max = 30 BR hybrid (in units of uwmap) is shown in the rightmost column. 




combined WMAP, Acbar and BOOMERanG analysis in 
terms of marginalized probability distributions. With a 
few exceptions, the low-£ bias results in a shift of a few 
tenths of a sigma. Perhaps most interesting is the effect 
on the much debated n a : This distribution is shifted to- 
wards higher values by the new low-£ likelihood, and the 
evidence for n s =/= 1 is thus reduced. 

Finally, we note that with Planck's reach to higher £'s, 
the 2 < £ < 30 data will not be critically important for 



constraining n s ; removing £ < 30 will increase a(n s ) by 
less than 10%. Planck's constraints on n s will thus not 
be dependent on low-£ polarization measurements, and 
therefore less sensitive to residual foregrounds. 

To summarize, we find that the improvements we 
make to the low-£ likelihood have a small but notice- 
ab le impact on the co smological parameters reported 
bv ISnererel et alJ ((2006). In particular, the evidence for 
n s 1 is weakened by 0.4a, which is significant given 
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the initial marginal ~ 2.5 — 3(7 detection. For full de- 
tails on results correcting for both this \ow-£ statisti- 
cal issue and the high-l point source issue, we r efer t o 
iHuffenberger et all <|2006[) and lKristiansen et all (|2006fl . 

7. CONCLUSIONS 

In this paper we have made an extensive re-analysis 
of the WMAP temperature anisotropy data. We have 
re-computed the power spectrum using six different 
codes, and the cosmological parameters with two differ- 
ent codes. Five different groups have contributed with 
separate computations. In total, the amount of cross- 
checks provided by this effort leads us to conclude that 
most important analysis issues concerning the three-year 
temperature power spectrum are well understood. 

Our main conclusion from this work is that we confirm 
the WMAP temperature power spectrum on most scales. 
However, subtle discrepancies are found both at large and 
small scales, and these can be summarized as follows: 

1. There is a small but significant bias at large angu- 
lar scales (5-10% at £ < 30) in the WMAP power 
spectrum and published likelihood code. This is 
primarily due to statistical issues in the form of 
an inadequate likelihood approximation, and sec- 
ondarily due to residual foregrounds. 

2. There is a systematic bias between the V- and 
W-bands. This h a s rece ntly been identified by 
Huff enberger et al . ( 2006) to be at least par- 
tially caused by over-correction of unresolved point 
sources. However, even after correcting for this, a 
small discrepancy is present, and this could possi- 
bly indicate further systematics. A possible candi- 
date worth considering could be uncorrected beam 
asymmetries. 

3. The low-^? likelihood bias has a noticable effect of 
some parameters. Most interestingly, it increases 
n s by 0.4ct, reducing the nominal detection of n s ^ 
1 from ~ 2.7 to ~ 2.3(7. Further, as reported by 
Huff enberger et aD l|2006j) . correcting for the point 
source over-correction increases n s by an additional 
~ 0.3cr, to a final significance of n s ^ 1 of ~ 2.0. 

In addition to these cosmologically important results, 
we make a few algorithmic recommendations based on 
the work presented here regarding the WMAP analysis. 
First and foremost, an exact likelihood evaluation should 
be used at least up to i < 30. Either direct evaluation 
or Gibbs sampling can be used for these scales, but of 
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course, the negligible computational cost of the latter, 
when incorporated into an MCMC sampler, makes it the 
preferred choice. Second, we advocate using the more 
conservative Kp2e cut at very low 1 as a hedge against 
possible foreground contamination just outside the di- 
rectly downgraded Kp2 mask. However, this has a neg- 
ligible effect in terms of cosmological parameters. 

The products from the Gibbs sampling analysis (with 
a uniform prior) are made publically available 14 in the 
form of three data sets. First, the most relevant for the 
applications presented in this paper is the collection of 
4 x 1000 sky signal spectrum samples. These form the ba- 
sis of the Blackwell-Rao estimator, and may be used for 
parameter estimation as demonstrated in this paper. A 
corresponding Fortran 90 module is also provided. Sec- 
ond, a set of 4000 ensemble averaged spectrum samples 
are presented. These may be used for visualization pur- 
poses of the power spectrum. Third, a set of sky map 
samples are provided for I < 50. These may be used 
for analyses that require phase information, for instance 
studies of non-Gaussianity. 

Thus, we conclude that although fast methods such as 
MASTER are very useful in the early stages of analyz- 
ing a new experiment, when fast turn around times are 
essential, the extensions to Gibbs sampling (as described 
in the Appendices) have now rendered an exact method 
quite tractable with currently available computing re- 
sources. Based upon our experience with these methods, 
we believe that Gibbs sampling can and will play a sig- 
nificant role in the future analysis of Planck data. 
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APPENDIX 

A. LOW-i LIKELIHOOD EVALUATION 

Low-l likelihood evaluation from high-resolution data has received considerable attention in the last few years 
ijEfstathiod 120041 ISlosar et alJl200l IHinshaw et alll200l. The reason is simply that while the currently popular 
pseudo-spectrum power spectrum estimators (|Hivonet^iLll2002|) work very well on intermediate and small angular 
scales, they are clearly sub-optimal on large scales. 

The log-likelihood for a pixelized Gaussian random field d, with vanishing mean and covariance matrix C, is given 

by 

-21og£ = d t C- 1 d + log|C|, (Al) 
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up to a normalization constant. Evaluation of this expression requires inversion of the pixel-pixel covariance matrix, 
and therefore scales as 0(Np ix ). Consequently, direct likelihood analysis of mega-pixel maps is not currently feasible. 

The covariance matrix consists of a series of terms depending on the data model, but at the very least one needs a 
signal term S, which for an assumed isotropic CMB is given by the angular power spectrum Ce, 

Sa = ^J2( 2i+ l)C^(cos%). (A2) 

1=2 

Here be is the Legendre transform of the (circularly symmetric) experimental beam, Pt(x) are the Legendre polynomials, 
9ij is the angle between pixels i and j, and l max is a sufficiently large multipole. The exact definition of "sufficiently 
large" will be made explicit shortly. 

Next, in most real-world applications there is also a term describing instrumental noise, N. Often, this is modeled 
as independent between pixels (white), and given by an overall noise level per observation <7o and a total number of 
observations per pixel N h s (i), such that Nij = <j$/ \J N h s {i)5ij. For current CMB experiments, such as WMAP, this 
is strongly sub-dominant to the signal term on large angular scales, and may in principle be omitted without significant 
loss of accuracy. However, this is not entirely true, since it can (and should) play an important role in regularizing the 
covariance matrix. 

Finally, it is u s eful to include a number of template terms over which amplitudes one wishes to marginalize 
(IBond et al]ll998t lsl osar et aTll2004|) . Consider for instance a given foreground template f whose spatial structure is 
known, but overall amplitude is unknown. Then the corresponding covariance matrix is given by the outer product 
F = f f*, and by assigning a sufficiently large uncertainty to this structure the corresponding mode is projected out 
from the data. In total, the covariance matrix may be written as 

C = S + N + ^A. ( F 4 , (A3) 

i 

where are numerically sufficiently large constants. 

Conceptually speaking, the above description completely defines the likelihood problem. However, for the particular 
problem of low-resolution analysis with high signal-to-noise data there is one particular issue that must be considered 
in greater detail in order to produce robust results, and that is the relationship between a finite pixelization, bandwidth 
limitation, and covariance matrix regularization. 

For the expansion in eauation lA2l to be valid, £ max must be chosen sufficiently large such that there is negligible power 
beyond this multipole. At the same time it must be smaller than the corresponding Nyquist frequency of the chosen 
pixelization. For HEALPix 15 the recommended upper limit on this value is ^ max = 2AT s id c (where N p i x = 127V s 2 ido ), 
but through various numerical techniques it can be pushed to £ max = 3 A^ide- Beyond this, one is almost assured to 
get nonsensical results, and one must therefore make sure that the experimental beam suppresses all signal beyond 
this value. (A good rule-of-thumb is that the beam FWHM must be at least 2.5 times the pixel size. For example, at 
-/Vside = 16 the pixel size is 3.6°, and the smallest beam supported by this pixelization is 9° FWHM.) 

In a pure signal map s, with a given bandwidth limit ^ max , there are (£ max + l) 2 real modes. Enforcing the 
recommended Nyquist limit of ^ max = 2 N a a c , a given pixelization can thus maximally resolve ~ 4A r s 2 ido independent 
modes, which is considerably less than the number of pixels in the map, Npi* = 12N^ ide . By a simple mode counting 
argument, it is clear that a signal-only covariance matrix necessarily will be singular, and consequently, it must be 
regularized in some way before inversion. 

The solution used by the currently available WMAP likelihood code is simply to increase £ max to 4 A^de, and thereby 
go beyond the Nyquist frequency of the pixelization. While it is true that this does indeed make the matrix invertible 
(since the number of projected modes is now greater than the number of pixels), it is also mathematically arbitrary, 
highly pixelization dependent, and not connected to the observed data. 

A much more reliable approach is readily available by means of the noise covariance matrix. By adding a small 
amount of white noise to the data, which has N p - lx independent modes, the matrix becomes invertible, and the data 
description is still accurate. A stable and well defined algorithm for computing low-resolution likelihoods from high- 
resolution data may be formulated as follows: 

1. Choose some £ max of interest. 

2. Determine the smallest pixelization that comfortably supports this scale and the corresponding pixel size 9 p - lx . 

3. Smooth the data with a Gaussian beam of 2.56* p ; x FWHM. 

4. Add Gaussian noise to the data with a variance such that the data are strongly noise dominated at 3N s id c - 

5. Compute the likelihood as described above, with appropriately defined covariance matrices. 



http:/ /healpix.jpl. nasa.gov/ 
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B. SOME GIBBS SAMPLING EXTENSIONS 

Pow er s pectrum estimation th rough Gibbs sampling was originally introduced to CMB applications b v I Jewell et al. 
(120041) andl Wandelt et afl J2004L and later implemented for high-resolution applications like WMAP bv lEriksenetaf 
(2004) and lOTJwvBr*eTld] ^2004) . The codes we use in the present paper are direct descendants from those described 
in the latter two papers, but with a few simple extensions we describe in this Appendix. Specifically, in order to reduce 
the Markov chain correlation length in the low signal-to-noise regime, we have implemented binned Cg sampling and 
sub-space sampling, and also sampling with Jeffrey's ignorance prior which has an effect on large angular scales. 

Let us first recall the basic idea of Gibbs sampling. Suppose we want to draw samples from a complicated joint 
probability distribution P(x,y), but only know how to sample from the conditional distributions P(x\y) and P(y\x). 
Then the theory of Gibbs sampling tells us that joint samples (x l ,y l ) may b e obtained by altern atel y drawing from 
these two distributions, x l+1 <— P(x\y l ) and y %+1 <— P(y\x t+1 ). As shown bv I Jewell et all l|2004j) and lWandelt et all 
(2004), this may be applied to CMB analysis because it is in fact feasible to sample from both P(Cg\d, s) and P(s|d, Cg) 
where d are observed data and s is the true CMB sky. Explicitly, the algorithm may be described by these two steps: 



s*+i«_P(s|C*,d), (Bl) 
C\ +1 ^P{Cg\s l+1 ). (B2) 

Here <— indicates sampling from the distribution on the right. Thus, after some burn-in time (Cj, s % ) will be drawn 
from the desired distribution, and these samples may subsequently be used to establish marginal P(C^|d) and P(s|d) 
if desired. 



Binned Cg sampling 

One of the major problems with the implementations described bv lEriksen et al.l (|2004f) was their inability to probe 
the low signal-to-noise regime. The reason was simply that the Markov chain step size between two consecutive Cg 
samples is given by cosmic variance alone, while the full posterior is given by both cosmic variance and noise. Thus, 
in order to move a significant distance in the high-i 7 regime, a very large number of steps is required. 

One way to improve on this situation is simply to bin the power spectrum, and thereby increase the signal-to-noisc 
ratio per sampled parameter. (Note that such binning must be done internally in the P(Cg |d,s) = P(Cg |s) sampler 
in order to earn an improvement - binning by post-processing does not have the same effect.) 

Before describing the binned sampling algorithm, it is useful to recall the single-£ algorithm for sampling from 
P(Cg\s): First compute the power spectrum of the signal map, and write it for convenience on the form o~£ = 
^2m=-e \ a em\ 2 - Next, draw 2^—1 Gaussian random variates ppg with zero mean and unit variance, and form the 
sum p\ = *Y^j=\ \Pi\ 2 - Tben the desired power spectrum sample is given by 

C? 1 = ^ (B3) 

Sampling binned Cgs is very similar. First, we define the binning scheme to be uniform in Cg£(£ +1), and choose 
some bin limits £ mm and £ max - We then form the quantity 

<7b= J2 ^+l)K™| 2 - (B4) 

The number of independent harmonic modes in this sum is n = (£ m &x + l) 2 — ^mm, an d therefore we draw n Gaussian 
random variates p 3 with zero mean and unit variance. Next, we form the sum 



The common bin sample value is then 



^ = £H 2 " (B5) 



Cb = — , (B6) 
Pb 



and the actual power spectrum sample coefficients are 

Cg = C h /£(£+!). (B7) 



Sub-space sampling 

A second technique to speed up convergence is sub-space sampling. As described above, Gibbs sampling simply 
means sampling one parameter at a time while conditioning on all others. If beneficial, we may therefore choose to 
sample only a few Cgs and cr^'s at a time while conditioning on all others. 
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Specifically, we may extend the basic sampling scheme given in Equation IB2I as follows. 



s iow P( s iow\Cl low , C} high , sjjigh, d), (B8) 
C^ow-m-Us^ 1 ), (B9) 

S high ^~ ^( S liigli|C],highi C!,low> sf^ 1 , d), (BIO) 



C^^P(C t ,w\4&). (Bll) 



Sampling from P(Ce\s) for a sub-set follows exactly the same algorithm as before. For P(s\ ov/ \Ce, Shigh, d) two trivial 
modifications must be made: The complementary sky signal that is conditioned upon must be subtracted from the 
data prior to sampling, and the corresponding Ce coefficients must be set to zero. 

The advantage of this partitioning lies in the relationship between pre-conditioning performance and correlation 
length: The Markov chain correlation length is very short in the high signal-to-noise regime but very long in the low 
signal-to-noise regime. Thus, in principle we would like to make a larger number of steps at high £'s than at low 
^'s, in order to obtain similar mixing properties everywhere. On the other hand, most of the computational expense 
for Gibbs sampling is spent on sampling from P(s|Q,d) for which a linear system must be solved using Conjugate 
Gradients. This linear system is dense in the high signal-to-noise regime, but nearly diagonal in the low signal-to-noise 
regime. Therefore, by conditioning on the computationally expensive high signal-to-noise components, we can sample 
the high-€ components more aggressively with a low computational cost per sample. 

For the analysis presented here, this is implemented through the following sampling scheme: 

• For each main sample stored on disk, 

— draw one all-scale sample (~ 120 CG iterations). 

— draw three intermediate-scale (£ > 300) samples (~ 20 CG iterations). 

* draw three small-scale (£ > 400) samples (~ 7 CG iterations) for each intermediate-scale sample. 

The pre-conditioner may be individually tuned to each case. In particular, an expensive pre-conditioner is used for 
the first case, and a cheap, diagonal pre-conditioner is used for the latter two. Also, if for a particular application one 
finds that convergence is slow for low £'s (such as foreground sampling), one ma y condition on , say, £ > 50 and solve 
the system exactly in one single iteration using the pre-conditioner described bv lEriksen et alJ l)2004|) . 

To summarize, it is straightforward to obtain good convergence even in the low signal-to-noise regime using Gibbs 
sampling with the introduction of binned Ce sampling and sub-space sampling. For the V-band analyses presented in 
this paper (iV s id e = 512, £ max = 700), we produced 4000 such de-correlated samples in three days with 16 processors, 
reaching a Gelman-Rubin statistic value of R — 1 < 0.05 for the last bin and R — 1 < 0.01 for £ < 600. 



Jeffrey's prior 

The Bayesian approaches to power spectrum estimation require an explicit statement of the prior probability dis- 
tribution of the power spectrum. This prior reflects the experimenter's knowledge, or lack thereof, about the power 
spectrum before the data are collected. 

The Jeffreys prior is often useful because it encodes a complete lack of knowledge about scale. If we have some 
parameter, such as a value of Ci, which we know is positive, but have no idea of its order of magnitude, then we can 
use P(Cg) — 1/Ce as our prior. This is the Jeffreys prior, and it gives equal weight to each logarithmic bin, reflecting 
the initial belief that Ce is equally likely to be in any of them. 

In practice, the scale of the parameter in question is not completely unknown. In the case of the CMB, for example, 
\fC~t < 2.7K, since the temperature of the CMB cannot be negative anywhere on the sky. While this information 
should technically be included in the prior, it is often not necessary to include it. In this case, the data already 
constrain the values of Ce so strongly that the above cutoff in the prior would have no effect on the final posterior 
distribution. 

The Jeffreys prior weights the posterior probability more toward low values of Ce than a uniform prior would. The 
uniform prior (P{Ce) = constant) is useful because the posterior is then exactly equal to the likelihood. We plot the 
posterior with both priors in figure 4, to show the effect of the Jeffreys prior in the final posterior distribution. 
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